﻿using System;

namespace MathService.Temp
{
	public static class Temp1
	{
		public static double F(double x) { return 1 / x; }

		public static void Test()
		{
			const double a = 0.5;
			const double b = 1.5;
			var M = Math.Log(b) - Math.Log(a);
			Console.WriteLine("число элементарных ортрезков?   ");
			var n = Convert.ToInt32(Console.ReadLine());

			var h = (b - a) / n;
			var s3 = .0;
			var s4 = .0;
			for(var i = 1; i <= n - 1; i++)
				if(i % 2 == 0)
					s4 += F(a + i * h);
				else
					s3 += F(a + i * h);
			var s6 = (s3 + s4);
			s3 = h * (s6 + F(a));
			s4 = h * (s6 + F(b));

			var oshibka = M - s3;
			Console.WriteLine("Метод  Интеграл       Ошибка");
			Console.WriteLine("");
			Console.WriteLine("прямоуг.лев.     {0,-10:F6} {1,11:E6}", s3, oshibka);
			oshibka = M - s4;
			Console.WriteLine("прямоуг.прав.    {0,-10:F6} {1,11:E6}", s4, oshibka);
			var s7 = h * (s6 + ((F(a) + F(b)) / 2));
			oshibka = M - s7;
			Console.WriteLine("трапеция         {0,-10:F6} {1,11:E6}", s7, oshibka);
			h = (b - a) / n;
			const double s1 = 0;
			double s2 = 0;
			for(var i = 1; i <= n - 1; i++)
				if(i % 2 == 0)
					s2 += F(a + i * h);
				else
					s3 += F(a + i * h);
			var s = (h / 3) * (F(a) + F(b) + (2 * s1) + (4 * s2));
			oshibka = M - s;
			Console.WriteLine("симпсон          {0,-10:F6} {1,11:E6}", s, oshibka);

			Console.ReadLine();
		}

		[Copyright("sasha_sinica", url = "http://www.cyberforum.ru/csharp-beginners/thread159307.html")]
		public static void ZeydelMethodTest()
		{
			const int n = 3;
			var A = new[,]
			{
				{-0.1,0.2,-0.1}, 
				{0.1,-0.2,-0.2}, 
				{-0.2,0.1,-0.1}
			};

			bool IsTrue;
			const double eps = 0.000001;
			var iter = 0;

			var X1 = new double[n];
			for(var i = 0; i < n; i++) X1[i] = 1;

			var X2 = new double[n];
			var B = new[] { 1.6, -2.3, 1.5 };
			do
			{
				for(var i = 0; i < n; i++)
				{
					double sum = 0;

					for(var j = 0; j < n; j++)
					{
						if(j != i)
							sum += A[i, j] * X1[j];

						X2[i] = (B[i] - sum) / A[i, i];

						X1[i] = X2[i];
					}
				}
				IsTrue = true;

				for(var i = 0; i < n; i++)
				{
					if(!(Math.Abs(X2[i] - X1[i]) > eps)) continue;
					IsTrue = false;
					break;
				}


				iter++;

			}
			while(!IsTrue);

			for(var i = 0; i < n; i++)
				Console.WriteLine("{0} \r\n", X2[i]);

			Console.WriteLine("{0} \r\n", iter);
		}

		public static void ZeydelMethodTest2()
		{
			double[] xo = { 4.1, 3.8 };
			var r = 0;
			var k = 0;
			var x = new double[2][].Initialize(i => new double[100]);
			const double E = 0.001;

			Func<double> dydx = () =>
			{
				if(r == 0)
					return 3 * x[0][k] * x[0][k] - 12 * x[1][k];
				if(r == 1)
					return 3 * x[1][k] * x[1][k] - 12 * x[0][k];
				throw new ApplicationException();
			};

			Func<int, int, double> dydxKR = (K, R) =>
			{
				if(R == 0)
					return 3 * x[0][K] * x[0][K] - 12 * x[1][K];
				if(R == 1)
					return 3 * x[1][K] * x[1][K] - 12 * x[0][K];
				throw new ApplicationException();
			};

			Func<double> d2ydx2 = () => 6 * x[r][k];

			Func<double> xNext = () => x[r][k] - dydx() / d2ydx2();
			x[0][0] = xo[0];
			x[1][0] = xo[1];

			for(k = 0; ; )
			{
				for(r = 0; r <= 1; )
				{
					if(Math.Abs(r == 1 ? dydxKR(k + 1, r) : dydx()) < E)
						r++;
					else
						x[r++][k + 1] = xNext();
				}

				var f = true;

				for(var i = 0; i <= 1; i++)
				{

					if(dydxKR(k + 1, i) > E)
						f = false;
				}

				if(f) break;
				k++;
			}


		}

		public static void ZeydelMethodtest3(Matrix A, double[] b, double[] x, double eps = 0.001)
		{
			var n = x.Length;
			if(b.Length != n)
				throw new ArgumentException("Некорректная длина вектора b. Не равна длине ветора x");
			if(A.N != n)
				throw new ArgumentException("Некорректная число строк матрицы A");



			// Условие окончания
			Func<double[], double[], bool> converge = (x_kurrent, x_last) =>
			{
				for(var i = 0; i < n; i++)
					if(Math.Abs(x_kurrent[i] - x_last[i]) >= eps)
						return false;
				return true;
			};

			/*
				Ход метода, где:
				a[n][n] - Матрица коэффициентов
				x[n], p[n] - Текущее и предыдущее решения
				b[n] - Столбец правых частей
				Все перечисленные массивы вещественные и
				должны быть определены в основной программе,
				также в массив x[n] следует поместить начальное
				приближение столбца решений (например, все нули)
			*/

			var X_current = x;
			var X_last = new double[n];

			do
			{
				for(var i = 0; i < n; i++)
					X_last[i] = X_current[i];

				for(var i = 0; i < n; i++)
				{
					double sum = 0;
					for(var j = 0; j < i; j++)
						sum += A[i, j] * X_current[j];
					for(var j = i + 1; j < n; j++)
						sum += A[i, j] * X_last[j];
					X_current[i] = (b[i] - sum) / A[i, i];
				}
			}
			while(!converge(X_current, X_last));
		}




		//< желаемая точность
		// N - размерность матрицы; A[N][N] - матрица коэффициентов, b[N] - столбец свободных членов,
		// X[N] - начальное приближение, ответ записывается также в X[N];
		public static void Jacobi(Matrix A, double[] b, double[] X, double eps = 0.001)
		{
			var N = X.Length;
			var x = new double[N];
			double norm; // норма, определяемая как наибольшая разность компонент столбца иксов соседних итераций.

			do
			{
				for(var i = 0; i < N; i++)
				{
					x[i] = b[i];
					for(var j = 0; j < N; j++)
						if(i != j) x[i] -= A[i, j] * X[j];
					x[i] /= A[i, i];
				}
				norm = Math.Abs(X[0] - x[0]);
				for(var i = 1; i < N; i++)
				{
					if(Math.Abs(X[i] - x[i]) > norm)
						norm = Math.Abs(X[i] - x[i]);
					X[i] = x[i];
				}
			} while(norm > eps);
		}
	}
}
